
//TESTING

int gauge_param(int spacetime_dim_size, int indx, int mu, int nc, int i, int j) 
{ 
	return ((((indx*spacetime_dim_size) + (mu % spacetime_dim_size))*nc + i)*nc + j);
}

__kernel void input_u(__global int *spacetime_dim_size, __global int *nc, __global float *m)
{
	int xidx = get_global_id(0);
	int mu = 0;
	for(mu = 0; mu < *spacetime_dim_size; mu++) {
		int i = 0;
		for(i = 0; i < *nc; i++) {
			int j = 0;
			for(j = 0; j < *nc; j++) {
				m[gauge_param(*spacetime_dim_size, xidx, mu, *nc, i, j)] = (i == j) ? 1.0f : 0.0f;
			}
		}
	}
}

float gaugeGet(int n, __global float *m, int idx, int nc, int mu, int i, int j) {
	if(mu>=0) return m[((idx*n + (mu % n))*nc+i)*nc+j];
	else return m[((idx*n + (-mu % n))*nc+j)*nc+i];
}

int pathSet(float c, __global float *m, int idx, int nc, int i, int j) {
	m[(idx*nc+i)*nc+j] = c;
	return 0;
}

float pathGet(__global float *m, int idx, int nc, int i, int j) {
	return m[(idx*nc+i)*nc+j];
}

int secondLoopPart1(int spacetime_dim_size, __global float *gauge_m, int idx, int nc, int mu, __global *path_m) {
	for(int i = 0; i < nc; i++) {
		for(int j = 0; j < nc; j++) {
			pathSet(gaugeGet(spacetime_dim_size, gauge_m, idx, nc, mu, i, j), path_m, idx, nc, i, j);
		}
	}
	return 0;
}

int secondLoopPart2(int spacetime_dim_size, int path_size, __global int *path, int idx, int nc, __global float *tmp, __global float *path_m, __global float *gauge_m) {
	for(int p=1; p<path_size; p++) {
		for(int i=0; i<nc; i++) {
			for(int j=0; j<nc; j++) {
				tmp[(i * nc) + j]=0;

				for(int k=0; k<nc; k++) {
					tmp[(i * nc) + j]+=pathGet(path_m, idx, nc, i, k)*gaugeGet(spacetime_dim_size, gauge_m, idx, nc, path[p], k, j);
				}

				pathSet(tmp[(i * nc) + j], path_m, idx, nc, i, j);
			}
		}
	}
	return 0;
}

__kernel void loop2(__global int *spacetime_dim_size, __global float *gauge_m, __global int *path_size, __global int *path, __global float *path_m, __global int *nc, __global float *tmp)
{
	int xidx = get_global_id(0);
	secondLoopPart1(*spacetime_dim_size, gauge_m, xidx, *nc, path[0], path_m);

	secondLoopPart2(*spacetime_dim_size, *path_size, path, xidx, *nc, tmp, path_m, gauge_m);
}